import warnings
warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import KFold, LeaveOneOut, cross_val_score, cross_val_predict
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.linear_model import Ridge, Lasso
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.preprocessing import StandardScaler
from scipy import ndimage
from skimage import filters, measure, morphology
import os

sns.set_style('whitegrid')
plt.rcParams['figure.dpi'] = 100

# Load data
images = np.load('images_data/images.npy')
targets = np.load('images_data/targets.npy')

print(f"Images: {images.shape}, dtype={images.dtype}")
print(f"Targets: {targets.shape}, dtype={targets.dtype}")
print(f"Target range: [{targets.min()}, {targets.max()}], mean={targets.mean():.1f}, std={targets.std():.1f}")
print(f"Unique target values: {len(np.unique(targets))}")

Images: (100, 1080, 1080, 6), dtype=uint8
Targets: (100,), dtype=int64
Target range: [3, 75], mean=42.7, std=11.1
Unique target values: 41

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

axes[0].hist(targets, bins=30, edgecolor='black', alpha=0.7, color='steelblue')
axes[0].set_xlabel('Target Value')
axes[0].set_ylabel('Count')
axes[0].set_title('Target Distribution')
axes[0].axvline(x=targets.mean(), color='red', linestyle='--', label=f'Mean={targets.mean():.1f}')
axes[0].legend()

axes[1].scatter(range(len(targets)), np.sort(targets), c='steelblue', s=20, alpha=0.7)
axes[1].set_xlabel('Sample Index (sorted)')
axes[1].set_ylabel('Target Value')
axes[1].set_title('Sorted Target Values')

plt.tight_layout()
plt.show()

# Visualize sample images across all 6 channels
channel_names = ['Ch0 (DNA?)', 'Ch1 (ER?)', 'Ch2 (RNA?)', 'Ch3 (AGP?)', 'Ch4 (Mito?)', 'Ch5 (BF?)']

# Pick 4 images with different target values
sample_idx = [np.argmin(targets), np.argmax(targets),
              np.argsort(targets)[25], np.argsort(targets)[75]]

fig, axes = plt.subplots(4, 6, figsize=(24, 16))

for row, idx in enumerate(sample_idx):
    for ch in range(6):
        img = images[idx, :, :, ch]
        vmax = np.percentile(img, 99.5)
        axes[row, ch].imshow(img, cmap='gray', vmin=0, vmax=max(vmax, 1))
        axes[row, ch].set_title(f'{channel_names[ch]}\nTarget={targets[idx]}', fontsize=9)
        axes[row, ch].axis('off')

plt.suptitle('Sample Images Across 6 Channels (4 different target values)', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()

# Per-channel statistics
print("Per-channel statistics:")
print(f"{'Channel':<12} {'Min':>6} {'Max':>6} {'Mean':>8} {'Std':>8} {'99th%':>8}")
print("-" * 52)
for ch in range(6):
    ch_data = images[:, :, :, ch]
    print(f"{channel_names[ch]:<12} {ch_data.min():>6} {ch_data.max():>6} "
          f"{ch_data.mean():>8.2f} {ch_data.std():>8.2f} {np.percentile(ch_data, 99):>8.1f}")

# Composite view of one image
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
idx = np.argsort(targets)[50]  # median target

# RGB composite of first 3 channels
rgb = np.stack([images[idx,:,:,0], images[idx,:,:,1], images[idx,:,:,2]], axis=-1)
rgb_scaled = (rgb.astype(float) / np.percentile(rgb, 99.5) * 255).clip(0, 255).astype(np.uint8)
axes[0].imshow(rgb_scaled)
axes[0].set_title(f'RGB Composite (Ch 0,1,2), Target={targets[idx]}')
axes[0].axis('off')

# Second RGB composite
rgb2 = np.stack([images[idx,:,:,3], images[idx,:,:,4], images[idx,:,:,5]], axis=-1)
rgb2_scaled = (rgb2.astype(float) / max(np.percentile(rgb2, 99.5), 1) * 255).clip(0, 255).astype(np.uint8)
axes[1].imshow(rgb2_scaled)
axes[1].set_title(f'RGB Composite (Ch 3,4,5), Target={targets[idx]}')
axes[1].axis('off')

# Channel 0 (likely DNA/nuclei) zoomed
axes[2].imshow(images[idx, 300:700, 300:700, 0], cmap='hot')
axes[2].set_title(f'Ch0 (Nuclei?) Zoomed, Target={targets[idx]}')
axes[2].axis('off')

plt.tight_layout()
plt.show()

Per-channel statistics:
Channel         Min    Max     Mean      Std    99th%
----------------------------------------------------

Ch0 (DNA?)        0    203     2.14     4.80     25.0

Ch1 (ER?)         0    255     7.43     6.54     34.0

Ch2 (RNA?)        0    147     7.66     5.19     27.0

Ch3 (AGP?)        0    117     4.09     3.43     18.0

Ch4 (Mito?)       0    255     5.87     6.45     33.0

Ch5 (BF?)         0    119    13.04     2.15     18.0

# Quick check: does mean intensity per channel correlate with target?
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

for ch in range(6):
    ax = axes[ch // 3, ch % 3]
    mean_intensities = images[:, :, :, ch].reshape(100, -1).mean(axis=1)
    ax.scatter(mean_intensities, targets, alpha=0.6, s=30, c='steelblue')
    corr = np.corrcoef(mean_intensities, targets)[0, 1]
    ax.set_title(f'{channel_names[ch]} (r={corr:.3f})')
    ax.set_xlabel('Mean Intensity')
    ax.set_ylabel('Target')

plt.suptitle('Mean Channel Intensity vs Target', fontsize=14)
plt.tight_layout()
plt.show()

def extract_features(img):
    """Extract features from a single 6-channel image."""
    features = {}
    
    for ch in range(6):
        channel = img[:, :, ch].astype(np.float64)
        prefix = f'ch{ch}'
        
        # Intensity statistics
        features[f'{prefix}_mean'] = channel.mean()
        features[f'{prefix}_std'] = channel.std()
        features[f'{prefix}_max'] = channel.max()
        features[f'{prefix}_p95'] = np.percentile(channel, 95)
        features[f'{prefix}_p99'] = np.percentile(channel, 99)
        features[f'{prefix}_skew'] = float(pd.Series(channel.ravel()).skew())
        
        # Fraction of bright pixels
        threshold = np.percentile(channel, 95)
        if threshold > 0:
            features[f'{prefix}_bright_frac'] = (channel > threshold).mean()
        else:
            features[f'{prefix}_bright_frac'] = 0.0
        
        # Entropy (texture)
        hist, _ = np.histogram(channel, bins=64, range=(0, 255))
        hist = hist / hist.sum()
        hist = hist[hist > 0]
        features[f'{prefix}_entropy'] = -np.sum(hist * np.log2(hist))
        
        # Gradient magnitude (edge content)
        gy, gx = np.gradient(channel)
        grad_mag = np.sqrt(gx**2 + gy**2)
        features[f'{prefix}_grad_mean'] = grad_mag.mean()
        features[f'{prefix}_grad_std'] = grad_mag.std()
    
    # Object counting on Channel 0 (likely DNA/nuclei)
    ch0 = img[:, :, 0].astype(np.float64)
    
    # Adaptive thresholding approach
    for t_name, threshold_val in [('otsu', None), ('p97', np.percentile(ch0, 97)), ('p95', np.percentile(ch0, 95))]:
        if t_name == 'otsu':
            try:
                threshold_val = filters.threshold_otsu(ch0)
            except ValueError:
                threshold_val = ch0.mean() + 2 * ch0.std()
        
        binary = ch0 > threshold_val
        # Clean up
        binary = morphology.remove_small_objects(binary, min_size=50)
        binary = morphology.binary_dilation(binary, morphology.disk(1))
        labeled = measure.label(binary)
        props = measure.regionprops(labeled)
        
        features[f'ch0_obj_count_{t_name}'] = len(props)
        if props:
            areas = [p.area for p in props]
            features[f'ch0_obj_mean_area_{t_name}'] = np.mean(areas)
            features[f'ch0_obj_total_area_{t_name}'] = np.sum(areas)
        else:
            features[f'ch0_obj_mean_area_{t_name}'] = 0
            features[f'ch0_obj_total_area_{t_name}'] = 0
    
    # Object counting on Channel 1 (ER/cytoplasm)
    ch1 = img[:, :, 1].astype(np.float64)
    try:
        t_otsu = filters.threshold_otsu(ch1)
    except ValueError:
        t_otsu = ch1.mean() + 2 * ch1.std()
    binary1 = ch1 > t_otsu
    binary1 = morphology.remove_small_objects(binary1, min_size=50)
    labeled1 = measure.label(binary1)
    features['ch1_obj_count'] = labeled1.max()
    
    # Cross-channel features
    for i in range(6):
        for j in range(i+1, 6):
            ci = img[:, :, i].astype(np.float64).ravel()
            cj = img[:, :, j].astype(np.float64).ravel()
            corr = np.corrcoef(ci, cj)[0, 1]
            features[f'corr_ch{i}_ch{j}'] = corr if not np.isnan(corr) else 0.0
    
    return features

# Extract features for all images
print("Extracting features from 100 images...")
all_features = []
for i in range(len(images)):
    if (i + 1) % 20 == 0:
        print(f"  Processing image {i+1}/100...")
    feats = extract_features(images[i])
    all_features.append(feats)

feature_df = pd.DataFrame(all_features)
print(f"\nExtracted {feature_df.shape[1]} features per image")
print(f"Feature matrix shape: {feature_df.shape}")
feature_df.head()

Extracting features from 100 images...

  Processing image 20/100...

  Processing image 40/100...

  Processing image 60/100...

  Processing image 80/100...

  Processing image 100/100...

Extracted 85 features per image
Feature matrix shape: (100, 85)

# Top correlated features with target
correlations = feature_df.corrwith(pd.Series(targets)).abs().sort_values(ascending=False)

print("Top 20 features most correlated with target:")
print("=" * 50)
for feat, corr in correlations.head(20).items():
    sign = '+' if feature_df[feat].corr(pd.Series(targets)) > 0 else '-'
    print(f"  {sign}{corr:.3f}  {feat}")

# Plot top 6 features vs target
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
top_feats = correlations.head(6).index.tolist()

for i, feat in enumerate(top_feats):
    ax = axes[i // 3, i % 3]
    ax.scatter(feature_df[feat], targets, alpha=0.6, s=30, c='steelblue')
    corr = feature_df[feat].corr(pd.Series(targets))
    ax.set_title(f'{feat}\n(r={corr:.3f})', fontsize=10)
    ax.set_xlabel('Feature Value')
    ax.set_ylabel('Target')
    # Add regression line
    z = np.polyfit(feature_df[feat], targets, 1)
    p = np.poly1d(z)
    x_line = np.linspace(feature_df[feat].min(), feature_df[feat].max(), 100)
    ax.plot(x_line, p(x_line), 'r--', alpha=0.7)

plt.suptitle('Top 6 Correlated Features vs Target', fontsize=14)
plt.tight_layout()
plt.show()

Top 20 features most correlated with target:
==================================================
  +0.947  ch0_obj_count_otsu
  +0.842  ch0_obj_count_p95
  +0.787  ch0_entropy
  +0.778  ch0_obj_total_area_otsu
  +0.767  ch0_mean
  -0.754  ch0_obj_mean_area_p95
  +0.753  ch0_obj_count_p97
  +0.744  ch0_p95
  +0.731  ch5_grad_mean
  -0.712  ch0_obj_mean_area_p97
  -0.701  ch4_skew
  +0.656  ch2_mean
  -0.618  ch0_skew
  +0.611  ch4_entropy
  +0.586  ch5_std
  +0.582  ch0_std
  +0.578  ch4_mean
  +0.560  ch3_entropy
  +0.547  ch0_p99
  +0.544  ch3_mean

# Prepare data
X = feature_df.values
y = targets

# Handle NaN/Inf
X = np.nan_to_num(X, nan=0.0, posinf=0.0, neginf=0.0)

# Scale features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Define models to evaluate
models = {
    'Ridge': Ridge(alpha=10.0),
    'Lasso': Lasso(alpha=1.0),
    'Random Forest': RandomForestRegressor(n_estimators=200, max_depth=5, random_state=42, n_jobs=-1),
    'Gradient Boosting': GradientBoostingRegressor(n_estimators=100, max_depth=3, learning_rate=0.1, random_state=42),
}

# 5-Fold CV (repeated 3x)
print("=" * 70)
print("5-FOLD CROSS-VALIDATION RESULTS (repeated 3x)")
print("=" * 70)
r2_hdr = "R^2"
print(f"{'Model':<25} {'MAE':>8} {'RMSE':>8} {r2_hdr:>8}")
print("-" * 55)

cv_results = {}

for name, model in models.items():
    all_mae, all_rmse, all_r2 = [], [], []
    
    for seed in [42, 123, 456]:
        kf = KFold(n_splits=5, shuffle=True, random_state=seed)
        
        y_pred_cv = cross_val_predict(model, X_scaled, y, cv=kf)
        mae = mean_absolute_error(y, y_pred_cv)
        rmse = np.sqrt(mean_squared_error(y, y_pred_cv))
        r2 = r2_score(y, y_pred_cv)
        
        all_mae.append(mae)
        all_rmse.append(rmse)
        all_r2.append(r2)
    
    mean_mae = np.mean(all_mae)
    mean_rmse = np.mean(all_rmse)
    mean_r2 = np.mean(all_r2)
    
    cv_results[name] = {'MAE': mean_mae, 'RMSE': mean_rmse, 'R2': mean_r2}
    print(f"{name:<25} {mean_mae:>8.2f} {mean_rmse:>8.2f} {mean_r2:>8.3f}")

# Identify best model
best_model_name = min(cv_results, key=lambda x: cv_results[x]['MAE'])
print(f"\nBest model by MAE: {best_model_name}")

======================================================================
5-FOLD CROSS-VALIDATION RESULTS (repeated 3x)
======================================================================
Model                          MAE     RMSE      R^2
-------------------------------------------------------

Ridge                         2.65     4.31    0.848
Lasso                         2.23     4.41    0.841

Random Forest                 2.54     5.07    0.789

Gradient Boosting             2.49     4.95    0.800

Best model by MAE: Lasso

# LOO CV on best model
print(f"Running Leave-One-Out CV for {best_model_name}...")
best_model = models[best_model_name]

loo = LeaveOneOut()
y_pred_loo = cross_val_predict(best_model, X_scaled, y, cv=loo)

mae_loo = mean_absolute_error(y, y_pred_loo)
rmse_loo = np.sqrt(mean_squared_error(y, y_pred_loo))
r2_loo = r2_score(y, y_pred_loo)

print(f"\nLOO CV Results ({best_model_name}):")
print(f"  MAE:  {mae_loo:.2f}")
print(f"  RMSE: {rmse_loo:.2f}")
print(f"  R^2:  {r2_loo:.3f}")

# Predicted vs Actual plot
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

axes[0].scatter(y, y_pred_loo, alpha=0.6, s=40, c='steelblue', edgecolors='white', linewidth=0.5)
min_val = min(y.min(), y_pred_loo.min()) - 2
max_val = max(y.max(), y_pred_loo.max()) + 2
axes[0].plot([min_val, max_val], [min_val, max_val], 'r--', alpha=0.7, label='Perfect prediction')
axes[0].set_xlabel('Actual Target')
axes[0].set_ylabel('Predicted Target')
axes[0].set_title(f'{best_model_name} - LOO CV\nMAE={mae_loo:.2f}, R^2={r2_loo:.3f}')
axes[0].legend()

# Residuals
residuals = y - y_pred_loo
axes[1].scatter(y, residuals, alpha=0.6, s=40, c='coral', edgecolors='white', linewidth=0.5)
axes[1].axhline(y=0, color='black', linestyle='-', alpha=0.3)
axes[1].set_xlabel('Actual Target')
axes[1].set_ylabel('Residual (Actual - Predicted)')
axes[1].set_title(f'Residual Plot\nMean residual={residuals.mean():.2f}')

plt.tight_layout()
plt.show()

Running Leave-One-Out CV for Lasso...

LOO CV Results (Lasso):
  MAE:  2.13
  RMSE: 4.24
  R^2:  0.853

# Train on full data for feature importance
if 'Random Forest' in models:
    rf = RandomForestRegressor(n_estimators=200, max_depth=5, random_state=42, n_jobs=-1)
    rf.fit(X_scaled, y)
    importances = rf.feature_importances_
    feat_names = feature_df.columns
    
    imp_df = pd.DataFrame({'feature': feat_names, 'importance': importances})
    imp_df = imp_df.sort_values('importance', ascending=False)
    
    fig, ax = plt.subplots(figsize=(12, 8))
    top20 = imp_df.head(20)
    ax.barh(top20['feature'][::-1], top20['importance'][::-1], color='steelblue')
    ax.set_xlabel('Feature Importance')
    ax.set_title('Top 20 Feature Importances (Random Forest)')
    plt.tight_layout()
    plt.show()
    
    print("\nTop 15 features:")
    for _, row in imp_df.head(15).iterrows():
        print(f"  {row['importance']:.4f}  {row['feature']}")

Top 15 features:
  0.5660  ch0_obj_count_otsu
  0.0900  ch0_obj_count_p95
  0.0540  ch0_skew
  0.0460  ch0_mean
  0.0410  ch0_obj_mean_area_p95
  0.0332  ch0_obj_total_area_otsu
  0.0159  ch4_skew
  0.0130  ch0_entropy
  0.0089  ch0_obj_total_area_p95
  0.0074  ch4_grad_mean
  0.0073  ch0_grad_mean
  0.0056  ch3_skew
  0.0049  corr_ch1_ch4
  0.0046  ch0_obj_mean_area_otsu
  0.0042  ch4_mean

# Direct nuclei counting using Channel 0 (DNA)
def count_nuclei(img_ch0, min_size=30, max_size=5000):
    """Count nuclei-like objects in a channel image."""
    ch = img_ch0.astype(np.float64)
    
    # Try multiple thresholds and return the best
    results = {}
    
    for p in [95, 96, 97, 98]:
        threshold = np.percentile(ch, p)
        if threshold < 1:
            continue
        binary = ch > threshold
        binary = morphology.remove_small_objects(binary, min_size=min_size)
        
        # Watershed-like separation: erode then dilate
        binary = morphology.binary_erosion(binary, morphology.disk(2))
        binary = morphology.binary_dilation(binary, morphology.disk(2))
        
        labeled = measure.label(binary)
        props = measure.regionprops(labeled)
        
        # Filter by size
        valid = [p for p in props if min_size < p.area < max_size]
        results[f'p{p}'] = len(valid)
    
    return results

# Count for all images
print("Counting nuclei across all images...")
count_results = []
for i in range(100):
    counts = count_nuclei(images[i, :, :, 0])
    counts['idx'] = i
    counts['target'] = targets[i]
    count_results.append(counts)

count_df = pd.DataFrame(count_results)

# Find best counting threshold
print("\nNuclei count correlations with target:")
for col in [c for c in count_df.columns if c.startswith('p')]:
    r = count_df[col].corr(count_df['target'])
    mae = mean_absolute_error(count_df['target'], count_df[col])
    print(f"  {col}: r={r:.3f}, MAE={mae:.1f}")

# Plot best counting method
best_count_col = max([c for c in count_df.columns if c.startswith('p')],
                     key=lambda c: abs(count_df[c].corr(count_df['target'])))

fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(count_df[best_count_col], count_df['target'], alpha=0.6, s=40, c='coral')
ax.plot([0, count_df[best_count_col].max()], [0, count_df[best_count_col].max()], 'r--', alpha=0.5)
r = count_df[best_count_col].corr(count_df['target'])
mae = mean_absolute_error(count_df['target'], count_df[best_count_col])
ax.set_xlabel(f'Counted Objects ({best_count_col} threshold)')
ax.set_ylabel('Target Value')
ax.set_title(f'Direct Object Counting vs Target\nr={r:.3f}, MAE={mae:.1f}')
plt.tight_layout()
plt.show()

Counting nuclei across all images...

Nuclei count correlations with target:
  p95: r=0.843, MAE=7.4
  p96: r=0.820, MAE=8.0
  p97: r=0.775, MAE=6.7
  p98: r=0.635, MAE=7.8

# Final summary
print("=" * 70)
print("FINAL RESULTS SUMMARY")
print("=" * 70)

print("\n1. DATA OVERVIEW:")
print(f"   - 100 images, 1080x1080, 6 channels (Cell Painting)")
print(f"   - Target range: [{targets.min()}, {targets.max()}], likely nuclei/cell counts")
print(f"   - Regression problem with small dataset")

print("\n2. EVALUATION STRATEGY:")
print(f"   - 5-Fold CV (3x repeated) for model comparison")
print(f"   - Leave-One-Out CV for final assessment")
print(f"   - Metrics: MAE, RMSE, R^2")

print("\n3. MODEL RESULTS (5-Fold CV):")
for name, res in sorted(cv_results.items(), key=lambda x: x[1]['MAE']):
    print(f"   {name:<25} MAE={res['MAE']:.2f}, RMSE={res['RMSE']:.2f}, R^2={res['R2']:.3f}")

print(f"\n4. BEST MODEL ({best_model_name}) LOO CV:")
print(f"   MAE={mae_loo:.2f}, RMSE={rmse_loo:.2f}, R^2={r2_loo:.3f}")

print("\n5. KEY FINDINGS:")
print(f"   - Top correlated features: {', '.join(correlations.head(3).index.tolist())}")
print(f"   - Object counting approach: r={r:.3f}")
print(f"   - Image-derived features can predict target with reasonable accuracy")

======================================================================
FINAL RESULTS SUMMARY
======================================================================

1. DATA OVERVIEW:
   - 100 images, 1080x1080, 6 channels (Cell Painting)
   - Target range: [3, 75], likely nuclei/cell counts
   - Regression problem with small dataset

2. EVALUATION STRATEGY:
   - 5-Fold CV (3x repeated) for model comparison
   - Leave-One-Out CV for final assessment
   - Metrics: MAE, RMSE, R^2

3. MODEL RESULTS (5-Fold CV):
   Lasso                     MAE=2.23, RMSE=4.41, R^2=0.841
   Gradient Boosting         MAE=2.49, RMSE=4.95, R^2=0.800
   Random Forest             MAE=2.54, RMSE=5.07, R^2=0.789
   Ridge                     MAE=2.65, RMSE=4.31, R^2=0.848

4. BEST MODEL (Lasso) LOO CV:
   MAE=2.13, RMSE=4.24, R^2=0.853

5. KEY FINDINGS:
   - Top correlated features: ch0_obj_count_otsu, ch0_obj_count_p95, ch0_entropy
   - Object counting approach: r=0.843
   - Image-derived features can predict target with reasonable accuracy

from sklearn.linear_model import HuberRegressor, TheilSenRegressor

robust_models = {
    'HuberRegressor (e=1.35)': HuberRegressor(epsilon=1.35, alpha=1.0, max_iter=500),
    'HuberRegressor (e=1.1)': HuberRegressor(epsilon=1.1, alpha=1.0, max_iter=500),
    'TheilSenRegressor': TheilSenRegressor(random_state=42, n_jobs=-1),
    'Lasso (baseline)': Lasso(alpha=1.0),
}

print("=" * 70)
print("STRATEGY 1: ROBUST REGRESSION - LOO CV")
print("=" * 70)

loo = LeaveOneOut()
robust_results = {}

for name, model in robust_models.items():
    y_pred = cross_val_predict(model, X_scaled, y, cv=loo)
    mae = mean_absolute_error(y, y_pred)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    r2 = r2_score(y, y_pred)
    robust_results[name] = {'MAE': mae, 'RMSE': rmse, 'R2': r2, 'y_pred': y_pred}
    print(f"  {name:<30} MAE={mae:.2f}, RMSE={rmse:.2f}, R2={r2:.3f}")

# Plot comparison
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
for ax, (name, res) in zip(axes.ravel(), robust_results.items()):
    residuals_r = y - res['y_pred']
    ax.scatter(y, residuals_r, alpha=0.6, s=40, edgecolors='white', linewidth=0.5)
    ax.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax.set_xlabel('Actual Target')
    ax.set_ylabel('Residual')
    ax.set_title(f'{name}\nMAE={res["MAE"]:.2f}, R2={res["R2"]:.3f}')
    # Add trend line to show bias
    z = np.polyfit(y, residuals_r, 1)
    p = np.poly1d(z)
    x_line = np.linspace(y.min(), y.max(), 100)
    ax.plot(x_line, p(x_line), 'r--', alpha=0.7, label=f'slope={z[0]:.3f}')
    ax.legend()

plt.suptitle('Strategy 1: Robust Regression - Residual Plots', fontsize=14)
plt.tight_layout()
plt.savefig('strategy1_robust_regression.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: strategy1_robust_regression.png")

======================================================================
STRATEGY 1: ROBUST REGRESSION - LOO CV
======================================================================

  HuberRegressor (e=1.35)        MAE=2.58, RMSE=4.32, R2=0.847

  HuberRegressor (e=1.1)         MAE=2.50, RMSE=4.37, R2=0.844

  TheilSenRegressor              MAE=5.99, RMSE=8.74, R2=0.376
  Lasso (baseline)               MAE=2.13, RMSE=4.24, R2=0.853

Saved: strategy1_robust_regression.png

from sklearn.compose import TransformedTargetRegressor

transformations = {
    'Lasso (no transform)': Lasso(alpha=1.0),
    'Lasso + log(1+y)': TransformedTargetRegressor(
        regressor=Lasso(alpha=1.0),
        func=np.log1p, inverse_func=np.expm1
    ),
    'Lasso + sqrt(y)': TransformedTargetRegressor(
        regressor=Lasso(alpha=1.0),
        func=np.sqrt, inverse_func=np.square
    ),
    'Ridge + log(1+y)': TransformedTargetRegressor(
        regressor=Ridge(alpha=10.0),
        func=np.log1p, inverse_func=np.expm1
    ),
}

print("=" * 70)
print("STRATEGY 2: TARGET TRANSFORMATION - LOO CV")
print("=" * 70)

transform_results = {}

for name, model in transformations.items():
    y_pred = cross_val_predict(model, X_scaled, y, cv=loo)
    mae = mean_absolute_error(y, y_pred)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    r2 = r2_score(y, y_pred)
    transform_results[name] = {'MAE': mae, 'RMSE': rmse, 'R2': r2, 'y_pred': y_pred}
    print(f"  {name:<30} MAE={mae:.2f}, RMSE={rmse:.2f}, R2={r2:.3f}")

# Plot comparison
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
for ax, (name, res) in zip(axes.ravel(), transform_results.items()):
    residuals_t = y - res['y_pred']
    ax.scatter(y, residuals_t, alpha=0.6, s=40, edgecolors='white', linewidth=0.5)
    ax.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax.set_xlabel('Actual Target')
    ax.set_ylabel('Residual')
    ax.set_title(f'{name}\nMAE={res["MAE"]:.2f}, R2={res["R2"]:.3f}')
    z = np.polyfit(y, residuals_t, 1)
    p = np.poly1d(z)
    x_line = np.linspace(y.min(), y.max(), 100)
    ax.plot(x_line, p(x_line), 'r--', alpha=0.7, label=f'slope={z[0]:.3f}')
    ax.legend()

plt.suptitle('Strategy 2: Target Transformation - Residual Plots', fontsize=14)
plt.tight_layout()
plt.savefig('strategy2_target_transform.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: strategy2_target_transform.png")

======================================================================
STRATEGY 2: TARGET TRANSFORMATION - LOO CV
======================================================================
  Lasso (no transform)           MAE=2.13, RMSE=4.24, R2=0.853
  Lasso + log(1+y)               MAE=8.54, RMSE=11.36, R2=-0.055
  Lasso + sqrt(y)                MAE=8.46, RMSE=11.22, R2=-0.028

  Ridge + log(1+y)               MAE=2.62, RMSE=3.79, R2=0.882

Saved: strategy2_target_transform.png

from sklearn.svm import SVR
from sklearn.model_selection import GridSearchCV

# --- SVR with hyperparameter tuning ---
print("Tuning SVR hyperparameters (5-Fold CV)...")
svr_param_grid = {
    'C': [1, 10, 50, 100],
    'epsilon': [0.1, 0.5, 1.0],
    'gamma': ['scale', 'auto']
}
svr_grid = GridSearchCV(SVR(kernel='rbf'), svr_param_grid, cv=5,
                        scoring='neg_mean_absolute_error', n_jobs=-1)
svr_grid.fit(X_scaled, y)
print(f"  Best SVR params: {svr_grid.best_params_}")
print(f"  Best 5-Fold MAE: {-svr_grid.best_score_:.2f}")

# --- XGBoost ---
try:
    from xgboost import XGBRegressor
    has_xgb = True
except ImportError:
    print("  XGBoost not installed, using GradientBoosting as substitute")
    has_xgb = False

nonlinear_models = {
    'SVR (tuned)': SVR(kernel='rbf', **svr_grid.best_params_),
    'SVR (linear)': SVR(kernel='linear', C=10),
    'Lasso (baseline)': Lasso(alpha=1.0),
}

if has_xgb:
    nonlinear_models['XGBoost'] = XGBRegressor(
        n_estimators=200, max_depth=3, learning_rate=0.05,
        subsample=0.8, colsample_bytree=0.8,
        reg_alpha=1.0, reg_lambda=1.0,
        random_state=42, verbosity=0
    )
else:
    nonlinear_models['GradientBoosting'] = GradientBoostingRegressor(
        n_estimators=200, max_depth=3, learning_rate=0.05,
        subsample=0.8, random_state=42
    )

print("\n" + "=" * 70)
print("STRATEGY 3: NON-LINEAR MODELS - LOO CV")
print("=" * 70)

nonlinear_results = {}

for name, model in nonlinear_models.items():
    y_pred = cross_val_predict(model, X_scaled, y, cv=loo)
    mae = mean_absolute_error(y, y_pred)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    r2 = r2_score(y, y_pred)
    nonlinear_results[name] = {'MAE': mae, 'RMSE': rmse, 'R2': r2, 'y_pred': y_pred}
    print(f"  {name:<30} MAE={mae:.2f}, RMSE={rmse:.2f}, R2={r2:.3f}")

# Predicted vs Actual
n_models = len(nonlinear_results)
fig, axes = plt.subplots(1, n_models, figsize=(6 * n_models, 5))
if n_models == 1:
    axes = [axes]

for ax, (name, res) in zip(axes, nonlinear_results.items()):
    ax.scatter(y, res['y_pred'], alpha=0.6, s=40, c='steelblue', edgecolors='white', linewidth=0.5)
    lims = [min(y.min(), res['y_pred'].min()) - 2, max(y.max(), res['y_pred'].max()) + 2]
    ax.plot(lims, lims, 'r--', alpha=0.7)
    ax.set_xlabel('Actual')
    ax.set_ylabel('Predicted')
    ax.set_title(f'{name}\nMAE={res["MAE"]:.2f}, R2={res["R2"]:.3f}')

plt.suptitle('Strategy 3: Non-linear Models - Predicted vs Actual', fontsize=14)
plt.tight_layout()
plt.savefig('strategy3_nonlinear_pred.png', dpi=150, bbox_inches='tight')
plt.show()

# Residual comparison
fig, axes = plt.subplots(1, n_models, figsize=(6 * n_models, 5))
if n_models == 1:
    axes = [axes]

for ax, (name, res) in zip(axes, nonlinear_results.items()):
    residuals_n = y - res['y_pred']
    ax.scatter(y, residuals_n, alpha=0.6, s=40, edgecolors='white', linewidth=0.5)
    ax.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax.set_xlabel('Actual Target')
    ax.set_ylabel('Residual')
    z = np.polyfit(y, residuals_n, 1)
    p = np.poly1d(z)
    x_line = np.linspace(y.min(), y.max(), 100)
    ax.plot(x_line, p(x_line), 'r--', alpha=0.7, label=f'slope={z[0]:.3f}')
    ax.set_title(f'{name}\nMAE={res["MAE"]:.2f}')
    ax.legend()

plt.suptitle('Strategy 3: Non-linear Models - Residual Plots', fontsize=14)
plt.tight_layout()
plt.savefig('strategy3_nonlinear_resid.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: strategy3_nonlinear_pred.png, strategy3_nonlinear_resid.png")

Tuning SVR hyperparameters (5-Fold CV)...
  Best SVR params: {'C': 50, 'epsilon': 0.1, 'gamma': 'auto'}
  Best 5-Fold MAE: 3.92

======================================================================
STRATEGY 3: NON-LINEAR MODELS - LOO CV
======================================================================
  SVR (tuned)                    MAE=3.76, RMSE=6.53, R2=0.651

  SVR (linear)                   MAE=3.11, RMSE=4.69, R2=0.820
  Lasso (baseline)               MAE=2.13, RMSE=4.24, R2=0.853

  XGBoost                        MAE=2.37, RMSE=4.50, R2=0.835

Saved: strategy3_nonlinear_pred.png, strategy3_nonlinear_resid.png

import torch
import torchvision.models as models
import torchvision.transforms as T

device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
print(f"Using device: {device}")

# Load pretrained ResNet50 as feature extractor (remove final FC layer)
resnet = models.resnet50(weights=models.ResNet50_Weights.IMAGENET1K_V1)
resnet = torch.nn.Sequential(*list(resnet.children())[:-1])  # Remove FC, keep avgpool
resnet = resnet.to(device)
resnet.eval()

# ImageNet normalization
normalize = T.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225])

def extract_resnet_features(imgs, channels_groups=[(0, 1, 2), (3, 4, 5)]):
    """Extract ResNet50 features from multi-channel images."""
    from PIL import Image
    all_features = []
    
    for i in range(len(imgs)):
        img_feats = []
        for group in channels_groups:
            rgb = np.stack([imgs[i, :, :, ch] for ch in group], axis=-1).astype(np.float32)
            for c in range(3):
                ch_max = rgb[:, :, c].max()
                if ch_max > 0:
                    rgb[:, :, c] = rgb[:, :, c] / ch_max
            
            pil_img = Image.fromarray((rgb * 255).astype(np.uint8))
            pil_img = pil_img.resize((224, 224), Image.BILINEAR)
            
            tensor_img = T.ToTensor()(pil_img)
            tensor_img = normalize(tensor_img)
            tensor_img = tensor_img.unsqueeze(0).to(device)
            
            with torch.no_grad():
                feat = resnet(tensor_img).squeeze().cpu().numpy()
            img_feats.append(feat)
        
        all_features.append(np.concatenate(img_feats))
        
        if (i + 1) % 20 == 0:
            print(f"  Processed {i+1}/{len(imgs)} images...")
    
    return np.array(all_features)

print("Extracting ResNet50 features from all 100 images...")
print("  (2 channel groups x 2048 dims = 4096 deep features per image)")
X_deep = extract_resnet_features(images)
print(f"  Deep feature matrix shape: {X_deep.shape}")

# Scale deep features
scaler_deep = StandardScaler()
X_deep_scaled = scaler_deep.fit_transform(X_deep)

Using device: cpu

Extracting ResNet50 features from all 100 images...
  (2 channel groups x 2048 dims = 4096 deep features per image)

  Processed 20/100 images...

  Processed 40/100 images...

  Processed 60/100 images...

  Processed 80/100 images...

  Processed 100/100 images...
  Deep feature matrix shape: (100, 4096)

# Evaluate deep features with different models
from sklearn.decomposition import PCA

# PCA to reduce 4096 dims (avoid curse of dimensionality with n=100)
for n_comp in [10, 20, 50]:
    pca = PCA(n_components=n_comp, random_state=42)
    X_pca = pca.fit_transform(X_deep_scaled)
    variance_explained = pca.explained_variance_ratio_.sum()
    print(f"PCA n={n_comp}: variance explained = {variance_explained:.3f}")

# Use PCA-reduced deep features
pca_best = PCA(n_components=20, random_state=42)
X_deep_pca = pca_best.fit_transform(X_deep_scaled)

# Also combine: hand-crafted + deep features
X_combined = np.hstack([X_scaled, X_deep_pca])
print(f"\nFeature sets:")
print(f"  Hand-crafted:  {X_scaled.shape}")
print(f"  Deep (PCA-20): {X_deep_pca.shape}")
print(f"  Combined:      {X_combined.shape}")

# Evaluate
deep_feature_sets = {
    'Lasso + HandCrafted (baseline)': (X_scaled, Lasso(alpha=1.0)),
    'Ridge + Deep (PCA-20)': (X_deep_pca, Ridge(alpha=10.0)),
    'Lasso + Deep (PCA-20)': (X_deep_pca, Lasso(alpha=1.0)),
    'Ridge + Combined': (X_combined, Ridge(alpha=10.0)),
}

print("\n" + "=" * 70)
print("STRATEGY 4: TRANSFER LEARNING FEATURES - LOO CV")
print("=" * 70)

deep_results = {}

for name, (X_feat, model) in deep_feature_sets.items():
    y_pred = cross_val_predict(model, X_feat, y, cv=loo)
    mae = mean_absolute_error(y, y_pred)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    r2 = r2_score(y, y_pred)
    deep_results[name] = {'MAE': mae, 'RMSE': rmse, 'R2': r2, 'y_pred': y_pred}
    print(f"  {name:<40} MAE={mae:.2f}, RMSE={rmse:.2f}, R2={r2:.3f}")

# Residual comparison
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
for ax, (name, res) in zip(axes.ravel(), deep_results.items()):
    residuals_d = y - res['y_pred']
    ax.scatter(y, residuals_d, alpha=0.6, s=40, edgecolors='white', linewidth=0.5)
    ax.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax.set_xlabel('Actual Target')
    ax.set_ylabel('Residual')
    ax.set_title(f'{name}\nMAE={res["MAE"]:.2f}, R2={res["R2"]:.3f}')
    z = np.polyfit(y, residuals_d, 1)
    p = np.poly1d(z)
    x_line = np.linspace(y.min(), y.max(), 100)
    ax.plot(x_line, p(x_line), 'r--', alpha=0.7, label=f'slope={z[0]:.3f}')
    ax.legend()

plt.suptitle('Strategy 4: Transfer Learning Features - Residual Plots', fontsize=14)
plt.tight_layout()
plt.savefig('strategy4_deep_features.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: strategy4_deep_features.png")

PCA n=10: variance explained = 0.491
PCA n=20: variance explained = 0.647
PCA n=50: variance explained = 0.865

Feature sets:
  Hand-crafted:  (100, 85)
  Deep (PCA-20): (100, 20)
  Combined:      (100, 105)

======================================================================
STRATEGY 4: TRANSFER LEARNING FEATURES - LOO CV
======================================================================
  Lasso + HandCrafted (baseline)           MAE=2.13, RMSE=4.24, R2=0.853
  Ridge + Deep (PCA-20)                    MAE=5.17, RMSE=6.88, R2=0.614
  Lasso + Deep (PCA-20)                    MAE=5.04, RMSE=6.51, R2=0.654

  Ridge + Combined                         MAE=2.66, RMSE=3.59, R2=0.894

Saved: strategy4_deep_features.png

# Collect all results into one table
all_results = {}

# Baseline
all_results['Lasso (baseline)'] = robust_results.get('Lasso (baseline)',
                                   transform_results.get('Lasso (no transform)'))

# Strategy 1: best robust
for name in ['HuberRegressor (e=1.35)', 'HuberRegressor (e=1.1)', 'TheilSenRegressor']:
    if name in robust_results:
        all_results[f'S1: {name}'] = robust_results[name]

# Strategy 2: best transform
for name in ['Lasso + log(1+y)', 'Lasso + sqrt(y)', 'Ridge + log(1+y)']:
    if name in transform_results:
        all_results[f'S2: {name}'] = transform_results[name]

# Strategy 3: best nonlinear
for name in ['SVR (tuned)', 'XGBoost', 'GradientBoosting']:
    if name in nonlinear_results:
        all_results[f'S3: {name}'] = nonlinear_results[name]

# Strategy 4: best deep
for name in ['Ridge + Deep (PCA-20)', 'Ridge + Combined']:
    if name in deep_results:
        all_results[f'S4: {name}'] = deep_results[name]

# Print sorted by MAE
print("=" * 75)
print("GRAND COMPARISON - ALL STRATEGIES (LOO CV)")
print("=" * 75)
print(f"{'Model':<45} {'MAE':>7} {'RMSE':>7} {'R2':>7}")
print("-" * 70)

sorted_results = sorted(all_results.items(), key=lambda x: x[1]['MAE'])
for name, res in sorted_results:
    marker = " *" if name == sorted_results[0][0] else ""
    print(f"  {name:<43} {res['MAE']:>7.2f} {res['RMSE']:>7.2f} {res['R2']:>7.3f}{marker}")

best_name = sorted_results[0][0]
print(f"\n* Best model: {best_name}")
print(f"  vs Lasso baseline: MAE improved by "
      f"{all_results['Lasso (baseline)']['MAE'] - sorted_results[0][1]['MAE']:.2f} "
      f"({(1 - sorted_results[0][1]['MAE']/all_results['Lasso (baseline)']['MAE'])*100:.1f}%)")

# Final visualization: best model vs baseline
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

baseline_pred = all_results['Lasso (baseline)']['y_pred']
best_pred = sorted_results[0][1]['y_pred']

for ax, (pred, title) in zip(axes, [(baseline_pred, 'Lasso (baseline)'),
                                      (best_pred, best_name)]):
    ax.scatter(y, pred, alpha=0.6, s=40, c='steelblue', edgecolors='white', linewidth=0.5)
    lims = [min(y.min(), pred.min()) - 2, max(y.max(), pred.max()) + 2]
    ax.plot(lims, lims, 'r--', alpha=0.7, label='Perfect')
    mae = mean_absolute_error(y, pred)
    r2 = r2_score(y, pred)
    ax.set_xlabel('Actual')
    ax.set_ylabel('Predicted')
    ax.set_title(f'{title}\nMAE={mae:.2f}, R2={r2:.3f}')
    ax.legend()

plt.suptitle('Baseline vs Best Model - Predicted vs Actual (LOO CV)', fontsize=14)
plt.tight_layout()
plt.savefig('grand_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: grand_comparison.png")

===========================================================================
GRAND COMPARISON - ALL STRATEGIES (LOO CV)
===========================================================================
Model                                             MAE    RMSE      R2
----------------------------------------------------------------------
  Lasso (baseline)                               2.13    4.24   0.853 *
  S3: XGBoost                                    2.37    4.50   0.835
  S1: HuberRegressor (e=1.1)                     2.50    4.37   0.844
  S1: HuberRegressor (e=1.35)                    2.58    4.32   0.847
  S2: Ridge + log(1+y)                           2.62    3.79   0.882
  S4: Ridge + Combined                           2.66    3.59   0.894
  S3: SVR (tuned)                                3.76    6.53   0.651
  S4: Ridge + Deep (PCA-20)                      5.17    6.88   0.614
  S1: TheilSenRegressor                          5.99    8.74   0.376
  S2: Lasso + sqrt(y)                            8.46   11.22  -0.028
  S2: Lasso + log(1+y)                           8.54   11.36  -0.055

* Best model: Lasso (baseline)
  vs Lasso baseline: MAE improved by 0.00 (0.0%)

Saved: grand_comparison.png

# Identify outliers from Lasso baseline
baseline_pred = robust_results['Lasso (baseline)']['y_pred']
residuals_baseline = y - baseline_pred

# Sort by absolute residual, pick top outliers
abs_resid = np.abs(residuals_baseline)
outlier_idx = np.argsort(abs_resid)[::-1]  # largest residual first

n_outliers = 5
print("Top outlier samples (largest |residual|):")
print(f"{'Rank':<6} {'Idx':<6} {'Actual':<10} {'Predicted':<12} {'Residual':<10}")
print("-" * 50)
for rank, idx in enumerate(outlier_idx[:n_outliers]):
    print(f"{rank+1:<6} {idx:<6} {y[idx]:<10} {baseline_pred[idx]:<12.1f} {residuals_baseline[idx]:<10.1f}")

# Also find a few "normal" samples for comparison
normal_idx = np.argsort(abs_resid)[:3]
print(f"\nBest-predicted samples (smallest |residual|) for comparison:")
for idx in normal_idx:
    print(f"  Idx={idx}, Actual={y[idx]}, Predicted={baseline_pred[idx]:.1f}, Residual={residuals_baseline[idx]:.1f}")

# --- Visualize outlier images (all 6 channels + RGB composite) ---
channel_names = ['Ch0 (DNA?)', 'Ch1 (ER?)', 'Ch2 (RNA?)', 'Ch3 (AGP?)', 'Ch4 (Mito?)', 'Ch5 (BF?)']

fig, axes = plt.subplots(n_outliers, 7, figsize=(28, 5 * n_outliers))

for row, idx in enumerate(outlier_idx[:n_outliers]):
    # 6 individual channels
    for ch in range(6):
        img_ch = images[idx, :, :, ch]
        vmax = np.percentile(img_ch, 99.5)
        axes[row, ch].imshow(img_ch, cmap='gray', vmin=0, vmax=max(vmax, 1))
        axes[row, ch].set_title(f'{channel_names[ch]}', fontsize=9)
        axes[row, ch].axis('off')
    
    # RGB composite (Ch0, Ch1, Ch2)
    rgb = np.stack([images[idx,:,:,0], images[idx,:,:,1], images[idx,:,:,2]], axis=-1)
    p995 = np.percentile(rgb, 99.5)
    if p995 > 0:
        rgb_scaled = (rgb.astype(float) / p995 * 255).clip(0, 255).astype(np.uint8)
    else:
        rgb_scaled = rgb.astype(np.uint8)
    axes[row, 6].imshow(rgb_scaled)
    axes[row, 6].set_title('RGB (Ch0-2)', fontsize=9)
    axes[row, 6].axis('off')
    
    # Row label
    axes[row, 0].set_ylabel(
        f'#{row+1} Idx={idx}\nActual={y[idx]}, Pred={baseline_pred[idx]:.1f}\nResid={residuals_baseline[idx]:.1f}',
        fontsize=11, fontweight='bold', rotation=0, labelpad=120, va='center'
    )

plt.suptitle(f'Top {n_outliers} Outlier Images (Largest |Residual| from Lasso Baseline)', fontsize=16, y=1.02)
plt.tight_layout()
plt.savefig('outlier_images.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: outlier_images.png")

# --- Now show "normal" well-predicted samples for contrast ---
fig2, axes2 = plt.subplots(len(normal_idx), 7, figsize=(28, 5 * len(normal_idx)))

for row, idx in enumerate(normal_idx):
    for ch in range(6):
        img_ch = images[idx, :, :, ch]
        vmax = np.percentile(img_ch, 99.5)
        axes2[row, ch].imshow(img_ch, cmap='gray', vmin=0, vmax=max(vmax, 1))
        axes2[row, ch].set_title(f'{channel_names[ch]}', fontsize=9)
        axes2[row, ch].axis('off')
    
    rgb = np.stack([images[idx,:,:,0], images[idx,:,:,1], images[idx,:,:,2]], axis=-1)
    p995 = np.percentile(rgb, 99.5)
    if p995 > 0:
        rgb_scaled = (rgb.astype(float) / p995 * 255).clip(0, 255).astype(np.uint8)
    else:
        rgb_scaled = rgb.astype(np.uint8)
    axes2[row, 6].imshow(rgb_scaled)
    axes2[row, 6].set_title('RGB (Ch0-2)', fontsize=9)
    axes2[row, 6].axis('off')
    
    axes2[row, 0].set_ylabel(
        f'Idx={idx}\nActual={y[idx]}, Pred={baseline_pred[idx]:.1f}\nResid={residuals_baseline[idx]:.1f}',
        fontsize=11, fontweight='bold', rotation=0, labelpad=120, va='center'
    )

plt.suptitle('Well-Predicted Samples (Smallest |Residual|) — For Comparison', fontsize=16, y=1.02)
plt.tight_layout()
plt.savefig('normal_images.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: normal_images.png")

# --- Residual plot with outliers highlighted ---
fig3, ax3 = plt.subplots(figsize=(10, 6))
ax3.scatter(y, residuals_baseline, alpha=0.5, s=40, c='steelblue', edgecolors='white', linewidth=0.5, label='Normal')
ax3.scatter(y[outlier_idx[:n_outliers]], residuals_baseline[outlier_idx[:n_outliers]], 
            s=120, c='red', edgecolors='black', linewidth=1, zorder=5, label='Top 5 outliers')
for rank, idx in enumerate(outlier_idx[:n_outliers]):
    ax3.annotate(f'#{rank+1} (idx={idx})', (y[idx], residuals_baseline[idx]),
                 textcoords="offset points", xytext=(8, 8), fontsize=9, color='red', fontweight='bold')
ax3.axhline(y=0, color='black', linestyle='-', alpha=0.3)
ax3.set_xlabel('Actual Target')
ax3.set_ylabel('Residual (Actual - Predicted)')
ax3.set_title('Residual Plot with Outliers Highlighted')
ax3.legend()
plt.tight_layout()
plt.savefig('residual_outliers_highlighted.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: residual_outliers_highlighted.png")

Top outlier samples (largest |residual|):
Rank   Idx    Actual     Predicted    Residual  
--------------------------------------------------
1      63     20         52.3         -32.3     
2      36     12         28.6         -16.6     
3      61     59         50.8         8.2       
4      93     64         56.9         7.1       
5      46     60         53.6         6.4       

Best-predicted samples (smallest |residual|) for comparison:
  Idx=74, Actual=33, Predicted=33.0, Residual=0.0
  Idx=39, Actual=44, Predicted=43.9, Residual=0.1
  Idx=18, Actual=33, Predicted=33.1, Residual=-0.1

Saved: outlier_images.png

Saved: normal_images.png

Saved: residual_outliers_highlighted.png

from sklearn.linear_model import HuberRegressor
from sklearn.compose import TransformedTargetRegressor

# Define candidate models
fix_models = {
    'Lasso (baseline)': Lasso(alpha=1.0),
    'Huber + sqrt(y)': TransformedTargetRegressor(
        regressor=HuberRegressor(epsilon=1.35, alpha=1.0, max_iter=500),
        func=np.sqrt, inverse_func=np.square
    ),
    'Huber + log(1+y)': TransformedTargetRegressor(
        regressor=HuberRegressor(epsilon=1.35, alpha=1.0, max_iter=500),
        func=np.log1p, inverse_func=np.expm1
    ),
    'Huber (no transform)': HuberRegressor(epsilon=1.35, alpha=1.0, max_iter=500),
    'Ridge + sqrt(y)': TransformedTargetRegressor(
        regressor=Ridge(alpha=10.0),
        func=np.sqrt, inverse_func=np.square
    ),
}

loo = LeaveOneOut()
fix_results = {}

print("=" * 70)
print("HUBER + SQRT ANALYSIS — LOO CV")
print("=" * 70)
print(f"  {'Model':<30} {'MAE':>7} {'RMSE':>7} {'R2':>7} {'Resid Slope':>12}")
print("-" * 70)

for name, model in fix_models.items():
    y_pred = cross_val_predict(model, X_scaled, y, cv=loo)
    mae = mean_absolute_error(y, y_pred)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    r2 = r2_score(y, y_pred)
    resid = y - y_pred
    slope = np.polyfit(y, resid, 1)[0]
    fix_results[name] = {'MAE': mae, 'RMSE': rmse, 'R2': r2, 'y_pred': y_pred, 'slope': slope}
    print(f"  {name:<30} {mae:>7.2f} {rmse:>7.2f} {r2:>7.3f} {slope:>12.4f}")

print(f"\n  (Slope closer to 0 = less systematic bias)")

# --- Residual comparison: 2x3 grid ---
fig, axes = plt.subplots(2, 3, figsize=(20, 12))
axes_flat = axes.ravel()

for i, (name, res) in enumerate(fix_results.items()):
    if i >= 6:
        break
    ax = axes_flat[i]
    resid = y - res['y_pred']
    
    ax.scatter(y, resid, alpha=0.6, s=40, edgecolors='white', linewidth=0.5,
               c='coral' if name == 'Lasso (baseline)' else 'steelblue')
    ax.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax.set_xlabel('Actual Target')
    ax.set_ylabel('Residual')
    
    # Trend line
    z = np.polyfit(y, resid, 1)
    p = np.poly1d(z)
    x_line = np.linspace(y.min(), y.max(), 100)
    ax.plot(x_line, p(x_line), 'r--', alpha=0.7, linewidth=2, label=f'slope={z[0]:.3f}')
    
    # Highlight outliers
    abs_r = np.abs(resid)
    top3 = np.argsort(abs_r)[-3:]
    ax.scatter(y[top3], resid[top3], s=100, c='red', edgecolors='black', linewidth=1, zorder=5)
    
    ax.set_title(f'{name}\nMAE={res["MAE"]:.2f}, R²={res["R2"]:.3f}, slope={res["slope"]:.3f}', fontsize=11)
    ax.legend(fontsize=10)

# Hide unused subplot if any
for j in range(len(fix_results), 6):
    axes_flat[j].set_visible(False)

plt.suptitle('Residual Comparison: Baseline vs HuberRegressor + Target Transform', fontsize=15, y=1.02)
plt.tight_layout()
plt.savefig('huber_sqrt_fix.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: huber_sqrt_fix.png")

# --- Predicted vs Actual: best fix vs baseline ---
best_fix_name = min(
    [(k, v) for k, v in fix_results.items() if k != 'Lasso (baseline)'],
    key=lambda x: abs(x[1]['slope'])
)[0]

fig2, axes2 = plt.subplots(1, 3, figsize=(21, 6))

for ax, name in zip(axes2, ['Lasso (baseline)', best_fix_name, 'Huber + sqrt(y)']):
    res = fix_results[name]
    ax.scatter(y, res['y_pred'], alpha=0.6, s=40, c='steelblue', edgecolors='white', linewidth=0.5)
    lims = [min(y.min(), res['y_pred'].min()) - 2, max(y.max(), res['y_pred'].max()) + 2]
    ax.plot(lims, lims, 'r--', alpha=0.7, label='Perfect')
    ax.set_xlabel('Actual')
    ax.set_ylabel('Predicted')
    ax.set_title(f'{name}\nMAE={res["MAE"]:.2f}, R²={res["R2"]:.3f}, slope={res["slope"]:.3f}')
    ax.legend()

plt.suptitle('Predicted vs Actual: Baseline → Fixed', fontsize=15)
plt.tight_layout()
plt.savefig('huber_sqrt_pred_vs_actual.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: huber_sqrt_pred_vs_actual.png")

======================================================================
HUBER + SQRT ANALYSIS — LOO CV
======================================================================
  Model                              MAE    RMSE      R2  Resid Slope
----------------------------------------------------------------------
  Lasso (baseline)                  2.13    4.24   0.853       0.2320

  Huber + sqrt(y)                   3.37    4.99   0.797       0.0606

  Huber + log(1+y)                  4.61    6.53   0.651       0.0929

  Huber (no transform)              2.58    4.32   0.847       0.0536
  Ridge + sqrt(y)                   2.43    3.81   0.881       0.1257

  (Slope closer to 0 = less systematic bias)

Saved: huber_sqrt_fix.png

Saved: huber_sqrt_pred_vs_actual.png

# ============================================================
# Step 1: Quantify Outlier Features
# ============================================================
# Get outlier indices from Lasso baseline
baseline_pred = robust_results['Lasso (baseline)']['y_pred']
residuals_baseline = y - baseline_pred
abs_resid = np.abs(residuals_baseline)
outlier_indices = np.argsort(abs_resid)[::-1][:5]

print("=" * 70)
print("OUTLIER FEATURE ANALYSIS")
print("=" * 70)

# Compute z-scores for all features
from scipy.stats import zscore
feature_z = pd.DataFrame(zscore(feature_df), columns=feature_df.columns)

for idx in outlier_indices[:3]:  # Top 3 outliers
    print(f"\n--- Sample Idx={idx} | Actual={y[idx]}, Predicted={baseline_pred[idx]:.1f}, Residual={residuals_baseline[idx]:.1f} ---")
    z_row = feature_z.iloc[idx]
    # Features where this sample is >2 std from mean
    extreme_feats = z_row[z_row.abs() > 2.0].sort_values(key=abs, ascending=False)
    if len(extreme_feats) > 0:
        print(f"  Features with |z-score| > 2.0 ({len(extreme_feats)} features):")
        for feat, zval in extreme_feats.head(10).items():
            actual_val = feature_df.iloc[idx][feat]
            pop_mean = feature_df[feat].mean()
            print(f"    {feat:<30} z={zval:>+6.2f}  (value={actual_val:.2f}, mean={pop_mean:.2f})")
    else:
        print("  No features with |z-score| > 2.0")

# Heatmap: top outliers vs population for key features
top_corr_feats = correlations.head(15).index.tolist()

fig, ax = plt.subplots(figsize=(16, 6))
compare_idx = list(outlier_indices[:5]) + list(np.argsort(abs_resid)[:3])  # outliers + normals
labels = [f'OUT idx={i}\ny={y[i]},r={residuals_baseline[i]:.0f}' for i in outlier_indices[:5]] + \
         [f'GOOD idx={i}\ny={y[i]},r={residuals_baseline[i]:.0f}' for i in np.argsort(abs_resid)[:3]]

data_for_heatmap = feature_z.iloc[compare_idx][top_corr_feats]
sns.heatmap(data_for_heatmap, annot=True, fmt='.1f', cmap='RdBu_r', center=0,
            xticklabels=top_corr_feats, yticklabels=labels, ax=ax,
            cbar_kws={'label': 'Z-score'})
ax.set_title('Feature Z-scores: Outliers vs Well-Predicted Samples\n(Top 15 correlated features)', fontsize=13)
plt.tight_layout()
plt.savefig('outlier_feature_heatmap.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: outlier_feature_heatmap.png")

======================================================================
OUTLIER FEATURE ANALYSIS
======================================================================

--- Sample Idx=63 | Actual=20, Predicted=52.3, Residual=-32.3 ---
  Features with |z-score| > 2.0 (9 features):
    ch0_obj_mean_area_otsu         z= +4.83  (value=3391.39, mean=2113.28)
    ch0_obj_total_area_otsu        z= +3.54  (value=172961.00, mean=90532.69)
    ch4_p99                        z= +2.85  (value=51.00, mean=31.47)
    corr_ch0_ch4                   z= +2.38  (value=0.33, mean=0.16)
    corr_ch2_ch4                   z= +2.35  (value=0.90, mean=0.66)
    corr_ch3_ch4                   z= +2.28  (value=0.93, mean=0.74)
    ch4_grad_std                   z= -2.20  (value=0.53, mean=1.31)
    corr_ch1_ch4                   z= +2.05  (value=0.88, mean=0.70)
    ch4_grad_mean                  z= -2.04  (value=0.33, mean=0.65)

--- Sample Idx=36 | Actual=12, Predicted=28.6, Residual=-16.6 ---
  Features with |z-score| > 2.0 (32 features):
    ch0_max                        z= +5.29  (value=203.00, mean=60.67)
    corr_ch0_ch4                   z= +4.98  (value=0.51, mean=0.16)
    ch1_max                        z= +4.29  (value=250.00, mean=73.06)
    ch1_grad_std                   z= +4.16  (value=1.05, mean=0.45)
    corr_ch0_ch5                   z= -4.11  (value=-0.13, mean=-0.01)
    corr_ch0_ch2                   z= +3.99  (value=0.56, mean=0.38)
    ch0_obj_mean_area_otsu         z= -3.77  (value=1117.38, mean=2113.28)
    ch0_skew                       z= +3.58  (value=11.54, mean=4.25)
    corr_ch0_ch1                   z= +3.43  (value=0.41, mean=0.18)
    ch0_grad_mean                  z= +3.42  (value=0.44, mean=0.23)

--- Sample Idx=61 | Actual=59, Predicted=50.8, Residual=8.2 ---
  Features with |z-score| > 2.0 (8 features):
    ch0_p99                        z= +3.01  (value=43.00, mean=22.42)
    ch0_std                        z= +2.95  (value=8.74, mean=4.52)
    ch0_mean                       z= +2.53  (value=4.03, mean=2.14)
    ch0_p95                        z= +2.53  (value=26.00, mean=13.27)
    ch0_grad_std                   z= +2.42  (value=1.04, mean=0.60)
    ch0_max                        z= +2.39  (value=125.00, mean=60.67)
    ch0_grad_mean                  z= +2.22  (value=0.37, mean=0.23)
    ch0_entropy                    z= +2.19  (value=1.09, mean=0.68)

Saved: outlier_feature_heatmap.png

# ============================================================
# Step 2: Systematic comparison of all optimization combos
# ============================================================
from sklearn.feature_selection import SelectKBest, f_regression
from sklearn.pipeline import Pipeline

loo = LeaveOneOut()

# --- C. Feature Selection: find optimal k ---
print("Finding optimal k for SelectKBest...")
k_scores = {}
for k in [10, 15, 20, 30, 50, 85]:
    if k > X_scaled.shape[1]:
        continue
    pipe = Pipeline([
        ('select', SelectKBest(f_regression, k=k)),
        ('model', Lasso(alpha=1.0))
    ])
    y_pred = cross_val_predict(pipe, X_scaled, y, cv=loo)
    mae = mean_absolute_error(y, y_pred)
    k_scores[k] = mae
    print(f"  k={k:<4} MAE={mae:.2f}")
best_k = min(k_scores, key=k_scores.get)
print(f"  Best k={best_k} (MAE={k_scores[best_k]:.2f})")

# --- Full comparison: all combos of {Lasso, Huber} x {raw, sqrt} x {all, selectk} ---
combos = {
    # Baseline
    'Lasso (original)': Pipeline([
        ('model', Lasso(alpha=1.0))
    ]),
    # A: sqrt transform only
    'Lasso + sqrt(y)': TransformedTargetRegressor(
        regressor=Lasso(alpha=1.0),
        func=np.sqrt, inverse_func=np.square
    ),
    # B: Huber only
    'Huber': Pipeline([
        ('model', HuberRegressor(epsilon=1.35, alpha=1.0, max_iter=500))
    ]),
    # C: Feature selection only
    f'Lasso + SelectK(k={best_k})': Pipeline([
        ('select', SelectKBest(f_regression, k=best_k)),
        ('model', Lasso(alpha=1.0))
    ]),
    # A+B: Huber + sqrt
    'Huber + sqrt(y)': TransformedTargetRegressor(
        regressor=HuberRegressor(epsilon=1.35, alpha=1.0, max_iter=500),
        func=np.sqrt, inverse_func=np.square
    ),
    # A+C: Lasso + sqrt + SelectK
    f'Lasso + sqrt + SelectK({best_k})': TransformedTargetRegressor(
        regressor=Pipeline([
            ('select', SelectKBest(f_regression, k=best_k)),
            ('model', Lasso(alpha=1.0))
        ]),
        func=np.sqrt, inverse_func=np.square
    ),
    # A+B+C: Huber + sqrt + SelectK (full combo)
    f'Huber + sqrt + SelectK({best_k})': TransformedTargetRegressor(
        regressor=Pipeline([
            ('select', SelectKBest(f_regression, k=best_k)),
            ('model', HuberRegressor(epsilon=1.35, alpha=1.0, max_iter=500))
        ]),
        func=np.sqrt, inverse_func=np.square
    ),
    # B+C: Huber + SelectK
    f'Huber + SelectK({best_k})': Pipeline([
        ('select', SelectKBest(f_regression, k=best_k)),
        ('model', HuberRegressor(epsilon=1.35, alpha=1.0, max_iter=500))
    ]),
}

print("\n" + "=" * 80)
print("FULL OPTIMIZATION COMPARISON — LOO CV")
print("=" * 80)
print(f"  {'Model':<38} {'MAE':>6} {'RMSE':>6} {'R2':>6} {'Slope':>7} {'Max|r|':>7}")
print("-" * 80)

combo_results = {}
for name, model in combos.items():
    y_pred = cross_val_predict(model, X_scaled, y, cv=loo)
    mae = mean_absolute_error(y, y_pred)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    r2 = r2_score(y, y_pred)
    resid = y - y_pred
    slope = np.polyfit(y, resid, 1)[0]
    max_resid = np.max(np.abs(resid))
    combo_results[name] = {
        'MAE': mae, 'RMSE': rmse, 'R2': r2,
        'y_pred': y_pred, 'slope': slope, 'max_resid': max_resid
    }
    print(f"  {name:<38} {mae:>6.2f} {rmse:>6.2f} {r2:>6.3f} {slope:>+7.3f} {max_resid:>7.1f}")

Finding optimal k for SelectKBest...
  k=10   MAE=2.06
  k=15   MAE=2.07
  k=20   MAE=2.09

  k=30   MAE=2.07
  k=50   MAE=2.10
  k=85   MAE=2.13
  Best k=10 (MAE=2.06)

================================================================================
FULL OPTIMIZATION COMPARISON — LOO CV
================================================================================
  Model                                     MAE   RMSE     R2   Slope  Max|r|
--------------------------------------------------------------------------------
  Lasso (original)                         2.13   4.24  0.853  +0.232    32.3
  Lasso + sqrt(y)                          8.46  11.22 -0.028  +1.011    39.4

  Huber                                    2.58   4.32  0.847  +0.054    27.8
  Lasso + SelectK(k=10)                    2.06   4.17  0.858  +0.222    32.3

  Huber + sqrt(y)                          3.37   4.99  0.797  +0.061    26.8
  Lasso + sqrt + SelectK(10)               8.46  11.22 -0.028  +1.011    39.4

  Huber + sqrt + SelectK(10)               1.86   3.99  0.870  +0.092    33.8

  Huber + SelectK(10)                      1.80   4.20  0.856  +0.112    36.0

# ============================================================
# Step 3: Head-to-head residual plot — Original vs Optimized
# ============================================================

# Pick key models to compare
show_models = ['Lasso (original)', 'Huber', 'Huber + sqrt(y)']
# Add best SelectK combo
selectk_models = [k for k in combo_results if 'SelectK' in k]
if selectk_models:
    best_sk = min(selectk_models, key=lambda k: abs(combo_results[k]['slope']))
    show_models.append(best_sk)

n_show = len(show_models)
fig, axes = plt.subplots(2, n_show, figsize=(7 * n_show, 12))

for col, name in enumerate(show_models):
    res = combo_results[name]
    pred = res['y_pred']
    resid = y - pred
    is_baseline = (name == 'Lasso (original)')
    color = 'coral' if is_baseline else 'steelblue'
    
    # Row 1: Predicted vs Actual
    ax1 = axes[0, col]
    ax1.scatter(y, pred, alpha=0.6, s=40, c=color, edgecolors='white', linewidth=0.5)
    lims = [min(y.min(), pred.min()) - 3, max(y.max(), pred.max()) + 3]
    ax1.plot(lims, lims, 'k--', alpha=0.4, label='Perfect')
    ax1.set_xlabel('Actual')
    ax1.set_ylabel('Predicted')
    ax1.set_title(f'{name}\nMAE={res["MAE"]:.2f}, R²={res["R2"]:.3f}', fontsize=11)
    ax1.legend()
    
    # Highlight outlier idx=63
    if 63 in range(len(y)):
        ax1.scatter(y[63], pred[63], s=150, c='red', edgecolors='black', linewidth=2, zorder=5, marker='*')
        ax1.annotate('idx=63', (y[63], pred[63]), textcoords="offset points",
                     xytext=(10, -10), fontsize=9, color='red', fontweight='bold')
    
    # Row 2: Residual plot
    ax2 = axes[1, col]
    ax2.scatter(y, resid, alpha=0.6, s=40, c=color, edgecolors='white', linewidth=0.5)
    ax2.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax2.set_xlabel('Actual Target')
    ax2.set_ylabel('Residual')
    
    # Trend line
    z = np.polyfit(y, resid, 1)
    p_line = np.poly1d(z)
    x_line = np.linspace(y.min(), y.max(), 100)
    ax2.plot(x_line, p_line(x_line), 'r-', alpha=0.8, linewidth=2.5, label=f'slope={z[0]:+.3f}')
    ax2.fill_between(x_line, p_line(x_line) - 0.5, p_line(x_line) + 0.5, alpha=0.1, color='red')
    
    # Highlight outliers
    top3 = np.argsort(np.abs(resid))[-3:]
    ax2.scatter(y[top3], resid[top3], s=100, c='red', edgecolors='black', linewidth=1.5, zorder=5)
    
    verdict = "BIASED" if abs(z[0]) > 0.15 else ("IMPROVED" if abs(z[0]) > 0.08 else "GOOD")
    ax2.set_title(f'Residuals: slope={z[0]:+.3f}, max|r|={res["max_resid"]:.1f}\n[{verdict}]',
                  fontsize=11, color='red' if verdict == "BIASED" else ('orange' if verdict == "IMPROVED" else 'green'))
    ax2.legend(fontsize=11)

plt.suptitle('Original vs Optimized: Predicted-vs-Actual (top) & Residuals (bottom)', fontsize=15, y=1.02)
plt.tight_layout()
plt.savefig('optimization_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: optimization_comparison.png")

# --- Summary verdict ---
print("\n" + "=" * 80)
print("OPTIMIZATION VERDICT")
print("=" * 80)
baseline = combo_results['Lasso (original)']
for name in show_models:
    res = combo_results[name]
    slope_change = abs(baseline['slope']) - abs(res['slope'])
    maxr_change = baseline['max_resid'] - res['max_resid']
    print(f"\n  {name}:")
    print(f"    Slope:     {baseline['slope']:+.3f} -> {res['slope']:+.3f}  (bias {'reduced' if slope_change > 0 else 'increased'} by {abs(slope_change):.3f})")
    print(f"    Max|resid|: {baseline['max_resid']:.1f} -> {res['max_resid']:.1f}  ({'better' if maxr_change > 0 else 'worse'} by {abs(maxr_change):.1f})")
    print(f"    MAE:       {baseline['MAE']:.2f} -> {res['MAE']:.2f}")
    print(f"    R²:        {baseline['R2']:.3f} -> {res['R2']:.3f}")

Saved: optimization_comparison.png

================================================================================
OPTIMIZATION VERDICT
================================================================================

  Lasso (original):
    Slope:     +0.232 -> +0.232  (bias increased by 0.000)
    Max|resid|: 32.3 -> 32.3  (worse by 0.0)
    MAE:       2.13 -> 2.13
    R²:        0.853 -> 0.853

  Huber:
    Slope:     +0.232 -> +0.054  (bias reduced by 0.178)
    Max|resid|: 32.3 -> 27.8  (better by 4.5)
    MAE:       2.13 -> 2.58
    R²:        0.853 -> 0.847

  Huber + sqrt(y):
    Slope:     +0.232 -> +0.061  (bias reduced by 0.171)
    Max|resid|: 32.3 -> 26.8  (better by 5.5)
    MAE:       2.13 -> 3.37
    R²:        0.853 -> 0.797

  Huber + sqrt + SelectK(10):
    Slope:     +0.232 -> +0.092  (bias reduced by 0.140)
    Max|resid|: 32.3 -> 33.8  (worse by 1.5)
    MAE:       2.13 -> 1.86
    R²:        0.853 -> 0.870

Metric	Value
Best Model	Huber Regressor + Sqrt Transform + SelectKBest
R-squared	$R^2 \approx 0.87$
MAE	$\approx 1.86$
Residual Quality	Randomly distributed around zero — no systematic bias

	ch0_mean	ch0_std	ch0_max	ch0_p95	ch0_p99	ch0_skew	ch0_bright_frac	ch0_entropy	ch0_grad_mean	ch0_grad_std	...	corr_ch1_ch2	corr_ch1_ch3	corr_ch1_ch4	corr_ch1_ch5	corr_ch2_ch3	corr_ch2_ch4	corr_ch2_ch5	corr_ch3_ch4	corr_ch3_ch5	corr_ch4_ch5
0	2.002874	4.255724	73.0	13.0	20.0	4.475237	0.049040	0.621754	0.191008	0.553978	...	0.738213	0.925405	0.603594	0.072857	0.822142	0.550737	0.103294	0.662061	0.084907	0.009577
1	2.211522	4.387540	75.0	13.0	22.0	4.429470	0.047583	0.647059	0.195780	0.617454	...	0.723932	0.937281	0.641773	-0.010399	0.855358	0.572129	0.045801	0.662521	0.017724	-0.050268
2	1.082278	3.427932	61.0	4.0	19.0	5.769195	0.047908	0.421443	0.173650	0.462299	...	0.775097	0.951925	0.682005	0.006831	0.837866	0.599758	-0.043919	0.715113	0.000344	-0.021342
3	1.534470	2.546892	37.0	8.0	13.0	3.966559	0.047555	0.519209	0.165943	0.352061	...	0.776422	0.923057	0.725164	-0.036261	0.877297	0.707432	-0.054379	0.765123	-0.049084	-0.088004
4	1.603420	2.911160	33.0	9.0	15.0	3.817244	0.048980	0.540433	0.183981	0.410624	...	0.828717	0.949342	0.773911	-0.093219	0.894399	0.724562	-0.134595	0.796943	-0.095269	-0.085401

Biological Image Analysis & Target Prediction¶

Analysis Goal¶

Challenges & Solutions¶

1. Systematic Prediction Bias¶

2. Outliers & Imaging Artifacts¶

3. Curse of Dimensionality & Overfitting¶

Key Results¶

Before & After: Residual Bias Correction¶

Part 1: Exploratory Data Analysis¶

1.1 Target Distribution¶

1.2 Image Visualization¶

1.3 Relationship Between Image Intensity and Targets¶

Part 2: Feature Engineering¶

2.1 Feature-Target Correlations¶

Part 3: Evaluation Strategy¶

Part 4: Model Development¶

4.1 Baseline Models¶

4.2 Leave-One-Out CV (Best Model)¶

4.3 Feature Importance¶

Part 5: Object Counting Approach¶

Part 6: Summary & Conclusions¶

Part 7: Addressing Systematic Bias — Four Improvement Strategies¶

Strategy 1: Robust Regression (HuberRegressor / TheilSenRegressor)¶

Strategy 2: Non-linear Target Transformation (log / sqrt)¶

Strategy 3: Non-linear Models (SVR / XGBoost)¶

Strategy 4: Transfer Learning Features (ResNet50)¶

Overall Comparison Across All Strategies¶

Outlier Analysis: Visualizing the Most Anomalous Samples¶

Recommended Fix: HuberRegressor + sqrt Target Transformation¶

Part 8: Systematic Optimization — Sqrt + Huber + Feature Selection¶

Step 1: Quantify Outlier Features¶

Step 2: Build Optimized Pipeline — Sqrt + Huber + SelectKBest¶

Step 3: Visual Diagnosis — Has the Bias Been Fixed?¶